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In inflationary cosmology, the particles constituting the Universe are created after inflation due 
to their interaction with moving inflaton field(s) in the process of preheating. In the fermionic 
sector, the leading channel is out-of equilibrium particle production in the non-perturbative regime 
of parametric excitation, which respects Pauli blocking but differs significantly from the perturbative 
expectation. We develop theory of fermionic preheating coupling to the inflaton, without and with 
expansion of the universe, for light and massive fermions, to calculate analytically the occupation 
number of created fermions, focusing on their spectra and time evolution. In the case of large 
resonant parameter q we extend for fermions the method of successive parabolic scattering, earlier 
developed for bosonic preheating. In an expanding universe parametric excitation of fermions is 
stochastic. Created fermions very quickly, within tens of inflaton oscillations, fill up a sphere of 
radius ~ g 1 ^ 4 in monetum space. We extend our formalism to the production of superheavy fermions 
and to 'instant' fermion creation. 
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I. INTRODUCTION 
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During cosmic inflation, it is assumed that entropy and temperature associated with particles of matter are diluted 
to practically zero values together with the number density of particles. After inflation, the inflaton field <j> oscillates 
around the minimum of its effective potential V(<p). The energy of inflaton oscillations is converted to the energy of 
newly created particles of matter. Eventually particles of different species are settled in a state of thermal equilibrium 
which marks the beginning of the conventional epoch of the hot Friedmann radiation domination. The actual process 
of particle creation from the classical background inflaton oscillations occurs very rapidly in the regime of parametric 
resonance (!]] before the thermal equilibrium will be settled. Particles are created non-perturbatively in the out- 
of-equilibrium state. The theory of this process, preheating, is elaborated in details for the creation of bosons [p]. 
For bosons (denoted x) > the leading effect is the stimulated process of particle creation in the regime of parametric 
• i-H , resonance, where the number density of created particles copiously increases with time as n x ~ exp J fidt. Soon, the 
backreaction of created x particles becomes important, so that the self-consistent dynamics of interacting bose fields 
H \ (f>{t,x) and x(t,x), which can be treated classically, can be revealed with lattice simulations [||. 

At the beginning of the preheating investigation, a study of fermion production did not look very interesting relative 
to the bosonic case. Indeed, the number density of fermions ip is bounded by Pauli blocking. Therefore, it was not 
expected that fermions will influence the dynamics of the inflaton and other scalar fields (despite a numerical study |J 
which down-played an observation that fermion production from an oscillating scalar is different from the perturbative 
prediction). 

However, it was understood with surprise |^^| that creation of fermions from the coherently oscillating inflaton field 
occurs very differently than what the conventional perturbation theory of inflaton decay <fr — ► ipip (say, due to a Yukawa 
coupling htjjcjyijj) would predict, and this may have many interesting cosmological applications. It turns out that the 
occupation number of fermions very quickly, within about ten inflaton oscillations, is saturated at the time-average 
value of about ~ 1/2. Moreover, in momentum space, fermions are excited within a non-degenerate "Fermi sphere" 
of a large radius k ~ q x l^m, where m is the frequency of inflaton oscillations, q is the usual dimensionless parameter 

of parametric excitation, q — — ^ . Ironically, it may be parametric excitations of fermions that are responsible for 
the most important observable signatures or observational constraints of preheating. Indeed, in some inflationary 
models there is significant production of gravitinos during the preheating stage |rj|j , Gravitinos are cosmologically 
dangerous relics, for a mass of ~ 100 Gev their abundance relative to that of the relic photons cannot exceed the 
bound n 3 / 2 / n ~( < 10~ 15 . Theory of gravitino production from preheating is rooted in the theory of spin 1/2 fermionic 
preheating. Another potentially important application of fermionic preheating is a possibility to produce superheavy 
fermions with a mass as large as ~ 10 18 Gev from inflatons of mass 10 13 Gev, as was noticed in Mm and investigated 
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in p] . Superheavy fermions may be interesting for the dark matter problem and for the problem of ultra-high energy 
cosmic rays. There are other interesting cosmological applications for the creation of fermions, e.g. the scenario of 
instant preheating IlH] and the creation of massive fermions during inflation Jlij ], Recently, fermionic preheating in 
hybrid inflation for some range of parameters was thoroughly studied |l2| . At present, it is hard to say how important 
fermionic preheating will be in the self-consistent non-linear dynamics of bose- and fermi-fields during preheating. 
Technical difficulty is to incorporate fermions into lattice simulations alongside with bosons. However, in (J] it was 
reported that even within the Hartree approximation backreaction of fermions catalyzes bosonic preheating. After 
all, for complete decay of inflatons one needs a 'three-legs' interaction, which is provided naturally by the Yukawa 
coupling with fermions. 

In this paper we develop the theory of fermionic preheating in an expanding universe, following our short paper ||, 
with an emphasis on the analytic results. In particular, we will generalize for the fermionic case some of the methods 
which we earlier developed for bosonic preheating [||^], in particular, the method of parabolic scattering, which works 
for large values of q, gives us an analytic formula for the occupation number of fermions nk(t) as function of time 
and momentum. We extend this method for production of superheavy fermions from a moving scalar field. We begin 
Sect. 2 of this paper with the Dirac equation and different VEVs for fermions interacting with a time-dependent 
background scalar field. In Sect. 3 the creation of fermions without expansion of the universe will be considered. We 
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mostly will consider scalar fields oscillating in a potential V — m ^ , although the methods can easily be applied to 
others. We will give a semi-analytic treatment of the problem based on some earlier results. For the case g> 1 we 
develop the method of successive parabolic scatterings. In Sect. 4 we take into account expansion of the universe, 
when preheating of fermions acquires new qualitative features. We extend the theory to describe the production of 
superheavy fermions. Our formalism also includes the case when fermions are created not from inflaton oscillations, 
but from a single instance, when inflaton field crosses a certain level. 



II. FORMALISM: FERMIONS COUPLING WITH BACKGROUND SCALAR IN FRW METRICS. 

We will consider the creation of spin-^ Dirac fermions i/j by a homogeneous, oscillating scalar field (j> in an expanding, 
flat FRW universe. Our strategy will be to solve the Heisenberg equation of motion for the quantum ?/>-fIeld in the 
presence of the classical backgrounds (f> and g^ v . Furthermore, we will assume the energy-momentum of the 0-field 
alone determines the expansion rate of the Universe. The matter action we use 



S M [<t>^,e^} = 



d A x e 



l <f> - V{4>) + iifrf D^ip- (rrty + htf>)iH> 



(1) 



contains a simple Yukawa coupling between the scalar field with effective potential V(<p) and fermion with bare- 
mass m^. Here 7 M is a space-time dependent Dirac gamma matrix, is the vierbein with e its determinant, and 
Df± = dfj, + ■g'JaP^u^ i s t ne spin-i covariant derrivative with vierbein dependent spin-connection, uo^ . The unbarred 
gammas are standard, Minkowski space-time Dirac matrices and "/ a p = JiaYp]- 



A. Classical Background Fields 

We consider flat FRW metrics ds 2 = dt 2 — a(i) 2 (ix 2 where a(t) is the scale factor of the universe. The Hubble 
parameter H = ~ is determined by the equation H 2 = 3 ^J a (^((f>) 2 + V^)^. The background homogeneous scalar 

field obeys the equation of motion c/) + 3H(j)+ ^ = 0. We will be most interested in fermion creation while the tfi- field 
oscillates about a minimum of its potential. For illustration of the methods, we will consider fermion creation in the 
context of chaotic inflation scenarios with potentials of the form V(4>) — jm 2 ^ 2 + -|</> 4 . After inflation, the 0-field 
oscillates quasi-periodically with a slowly decreasing amplitude due to the Hubble expansion. The specific value of the 
parameter or A (or some combination thereof) will be dictated by cosmology. We mostly will consider quadratic 
potential with A = 0. In this paper we ignore backreaction of created particles. 

B. Quantum Dirac Field 

Variation of the action ([!]) with respect to ■0 leads to the general relativistic generalization of the Dirac equation 

[vfDn - (m + h4f>{t))\ i>{x) = . (2) 
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It can be shown that for an FRW space-time we have 7 = 7 and 7* = jjjy7\ where the 7 Q 's are standard, Minkowski 
space-time Dirac matrices. Furthermore, spin connection leads to 



Thus, in our case, the Dirac equation becomes 



ij°d +i-i- V + 4 ) 7 - {m^ + hc/)) 



ip(x) = . 



(3) 



(4) 



Only standard gamma matrices now appear. Similarly, -ip — 1/^7° is found to obey the conjugate of Equation (Q). 

In the Heisenberg representation of QFT, the Dirac equation (||) or (ji|) becomes the equation of motion for the 
ip(x) field operator. ip(x) may be decomposed into eigen-spinors 



f (Pk / 

^,t) = J2 J j2nf ( a ^ Ufc -W e+Jk ' x + ^L v ^W' 



(•5) 



where u kl ±(t) is a positive-frequency eigenspinor of the Dirac equation (^) with helicity ±^ and v ki ±(t) = Cu|" ± (i) 
is its charge conjugate. Here C — 17072 is the standard charge conjugation matrix. To construct these eigen-spinors 
we use the ansatz 



Ufc,±(i) e 



+ik-x 



- i7 °a - i-7 • V - 4 ( - ) 7° - ("V + M 
a 2 V a ' 



Xfc(t)i2±(k)e 



+ik-x 



(G) 



where i?±(k) are eigenvectors of the helicity operator k • £ such that k • Si?±(k) = ±1 and 7°i?±(k) = +1. The 
time dependence of the eigen-spinor is contained intirely in the mode function Xk(t)- Substituting this ansatz into 
equation (Eh leads to the mode equation 



Xk + 4-Xfc 
a 



k 1 2 (aM eS ) 9 fa 



4 V a 



3 a 
2a 



Xk = 



where k 2 = |k| 2 and the effective mass of the fermions is 

M e ff = (m^j + hcj>) . 



(7) 



(8) 



The damping term in this equation may be removed by defining a new mode function X k (t) = a?Xk{t)- The mode 
equation becomes 



X k + 



k , ,1 (a M cS ) ■ A . ' 
a a 



X k = 



(9) 



where A(a) = (^) — \ (-)]• When a cx £™ such as in a matter or radiation dominated universe, A(a) ~ t~ 2 and 
can be neglected soon after inflation. 

Using the mode equation and ansatz, we find 



Ufc,±(£) = a 
Taking the charge conjugate of u k ,±, we find 

Vfc,±(t) = a -1 



iX k - 4 ( - I X k + (7 • k - M eff )X fc 



2 V a 



i?±(k) 



i?±(k) , 



(10) 



(11) 



where i?±(k) = —i7 2 i?:!j_(k) is an eigenvector of helicity such that k • Si?±(k) = ±1 and j°R±(k) = — 1. 

The energy-momentum tensor is obtained from the (i/> and tp symmetrized) matter action by variation with respect 
to the vierbein 



~~ 2 



(12) 
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and the Hamiltonian operator is 



H D = d 



.73 , 



(13) 



In general, if this Hamiltonian is diagonal in the annihilation (and creation) operators a k . s (at J and bk, s (&t s ) at 
t = 0, it will not be for later times. This is the signature of particle creation due to the time dependent background. 
In order to determine the number of particles produced, we perform a Bogliubov transformation on the creation and 
annihilation operators so as to diagonalize the Hamiltonian at time t. 

For the problem of particle creation quite often it is useful to represent the wave function in the adiabatic ( 
semi-classical, WKB) form 

X k (t) =a k N + e~ i fo dtnk(t) +f3 k N-e +% Jo * n *W > (14) 

where N± = (2f2fc(f2fc ± M e g)) -1 / 2 and the coefficients a k and (3 k correspond to the coefficients of the Bogliubov 
transformation. 

Once the Bogliubov transformation is done, we may write the comoving number density of particles rik(t) = |/3/c| 2 
in a given spin state through the solutions of the mode equation (^J) 

n k (t) = a ( Qk ~ M ° S ) \\X k \ 2 + Q 2 \X k \ 2 - 2Sl k Im{X k X* k )] , (15) 



2n k 

where M c g = (rn$ + htp) as before and il 2 = |y + Af 2 ff . The energy density in these particles is then 

2 f d 3 k If 
PH,{t) = t / -—r^nkit) = — / dkk 2 Q k n k (t) . (16) 



a 3 J (2tt) 2 

The normalization of the solutions X k (t) is such that X k (t — > Q _ ) = N+e~ inkt and n k (0) — 0. Thus, we find 
N + = (2Q k (£l k — M c ff)) -1 / 2 . These are the so-called positive frequency initial conditions. 

A comment about the regularization of fcrmion VEVs. In principle, fermionic VEVs like (T^ u ) require regularization 
, which in the presence of background metrics and scalars is rather non-trivial, see e.g. |13j and references therein. 
We will consider only the processes of particle creation, ignoring vacuum polarization. The creation of particles, 
which corresponds to the imaginary part of the effective action, has no formal divergencies and does not require 
regularization. Assuming the particle creation process dominates over vacuum polarization, we will not consider the 
issues of regularization. 

III. FERMIONIC PREHEATING WITHOUT EXPANSION OF THE UNIVERSE 

It is convenient to begin the investigation of fermionic preheating due to an oscillating scalar field with a simplified 
setting neglecting the expansion of the universe. This setting may have not only methodological advantages. Indeed, 
whenever the frequency of <f> oscillations is much greater than the rate of cosmic expansion H, it is sensible to neglect 
the time dependence of the scale factor in solving the mode equation (^|). This can be the case, for example, when 
spontaneous symmetry breaking occurs rapidly leaving the 0-field oscillating about a new minimum where the effective 
mass TO^ happens to be much larger than the Hubble parameter at the time of the transition. One such example is 
hybrid inflation scenarios which end with TeV scale energy densities. Another example is the inflationary model with 
the potential \<f) A . This theory possesses conformal properties: at the stage of inflaton oscillations equations for the 
fields by means of conformal transformations can be reduced to the equations in Minkowskii space-time, see e.g. M. 

In this paper we will mostly use a chaotic inflationary model with quadratic potential V{<f)) = ^m 2 ^ 2 . If we make 
the replacement a = 1 in all the formulas of Section 2, all the effects of expansion will be removed. Background 
oscillations take the form of harmonic oscillations <f>(t) — <fiof(t) with f(t) — cos(m^t) and </>o the time independent 
amplitude. It is convenient to define a new, dimensionless time variable r = m$t. With this change of variables, the 
mode equation (||) may be written 

X'l + [k 2 + (m +^/df) 2 - i^qf] X k = , (17) 

where we have introduced the dimensionless momentum k = the dimensionless fcrmion mass m= Hbt an d the 
resonance parameter q = — These three parameters completely determine the strength of the effect. In fact, this 
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form of the mode equation is valid not only for harmonic background oscillations, but for generic <f> oscillations in a 
general potential. For this, we identify with the frequency of oscillation, <p with its amplitude, and /(r) with the 
periodic background oscillations normalized to unit amplitude. Note that the frequency will be amplitude-dependent 
for a general, non-quadratic potential. 



A. Parametric Excitation of Fermions 



If individual inflatons at rest are decaying into light fermions in a process <f> ~~ > V'V'i perturbative calculations give 
the rate of decay ~ and the spectrum of created fermions is sharply peaked around m/2 with the width 

j, . In Figs. 0, |^, however, we plot the time-dependence and spectrum of occupation number for fermions created 
from the coherently oscillating inflaton field, as follows from numerical solution of Eq. (|17|) and ([L5J) . 
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FIG. 1. The occupation number rik of fermions in 
m 2 ^ 2 -theory as function of time ^ for a range of(q;n): 
dotted curve for (1.0; 0.05), bold curve for (10" 4 ;0.5), 
smooth sine-like curve for (10~ 2 ;0.5), utiggly sine-like 
curve for (1; 1.3). 
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FIG. 2. The comoving occupation number of fermions 
in m^if) 2 -inflation as a function of k 2 after 10 inflaton 
oscillations for resonance parameter q — 10. Also shown 
is the envelope function Fk defined by equation Jg"?[). 



The spectrum and time evolution are drastically different from what is expected from the perturbative calculations. 
We therefore can talk about a specific phenomena, the parametric excitation of fermions interacting with coherent 
background oscillations. 
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FIG. 3. Time dependence of the occupation num- 
ber in \cj> 4 -theory for q = 3 for two different modes: 
one is just inside the resonance band which grows like 
sinh(/ifct) ~ e Mfcr and another outside which oscillates 
like sin(/ifct) . 



FIG. 4. Bosonic occupation number after 10 back- 
ground oscillations for \(f> 4 -theory when q = = 3. A 
stripe with sharp edges corresponds to the e' lfeT instability. 



It is instructive to compare the spectrum and evolution of fcrmionic occupation number with those of bosonic 
occupation number. The mode function of a quantum bose scalar field x coupling as g 2 x' 2 4' 2 with oscillating inflaton 
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obeys the bosonic oscillator-like equation 



Y" 



2 



K + m b + q b f 



X bM = , (18) 



where q b = g m t" , m b = ■ In Figs. || and || we plot the spectrum and time evolution of bosonic occupation number 
calculated from bosonic oscillator-like equations (|l^). We use the familiar model of bosonic resonance due to the 
self-interaction in A0 4 inflation, which corresponds to g 2 — 3A, q b — 3, m b = 0, / (t) is given by oscillations in X(j) 
theory. 

In the bosonic case, there are distinct resonance bands in which bosonic modes are exponentially unstable, n b k (t) ~ 
gMfci" Particle creation takes place outside of the resonant bands as well, although there the occupation number is 
bounded and oscillates periodically. In the fermionic case rik < 1 is always bounded by Pauli blocking and oscillating 
periodically with time. The occupation number in both cases is changing in time. It is therefore convenient to 
introduce an envelope function Fk of the particle spectrum Q, which corresponds to averaged over short-time 
intervals (order of the background oscillation period) . A bosonic envelope function cannot be defined for the resonant 
bands (where it is e PfeT .) This zone corresponds to the gap between almost vertical lines in Fig. p. The fermionic 
envelope function Fk, on the contrary, can be defined everywhere. Although it is always bounded by unity, their 
structure is reminiscent of the resonant-band structure: for some range of k, the 'resonant' band, rik is close to unity, 
while in other ranges, the stable bands, it is significantly smaller if not zero. In jl2| different levels of Fk were plotted 
on the parameter plane (k, q), revealing a structure which reminiscent of the stability/instability chart of the bosonic 
parametric resonance. One of the most important results is that fermionic parametric excitations occurs very quickly, 
within about ten(s) background oscillations. Interestingly, fermionic "resonant bands" are excited the last, while 
non-resonant intervals fill first. 

Comparison of bosonic and fermionic parametric excitations is useful to understand some features of bosonic 
preheating. As we have seen, there is production of bosons outside of the resonance band. If we take into account the 
expansion of the universe, in the most interesting case of large q b , the difference between resonant and non-resonant 
excitations of bosons will be erased, and the regime of stochastic resonant production of bosons will be settled down 



B. Some Generic Analytic Results 



Dynamics of the Fermi field coupling to the background homogeneous scalar can be revealed with the second-order 
oscillator-like equation ( |l7| ) for the mode function Xf.(t). For periodic background oscillations f(t) = f(t + T), T is a 
period, some generic analytic results were derived a long time ago flo[| in the context of particle creation in a periodic 
external electromagnetic field. In particular, the occupation number of created particles at instances t — N S T, i.e. 
exactly after N s background oscillations, is given by expression jl0| 

n k {N s T) = ^^4^{lrnX^\T))\ (19) 
2,1 l k sin dk v ' 

where cosdk — ReXjp (T). To derive this result, one introduces two fundamental solutions of Eq. (|l7|), X^\t) and 

Xf \t); such that initially xj? 3 (0) = 1, a£ 1} (0) = and xf } (0) = 0, xf } (0) = 1. Expression @ involves only the 

value of the first fundamental solution XP^(T) exactly after the first oscillation. It says that the occupation number 
of created particles after N s background oscillations is modulated with a certain frequency Vk (which, as we will see, 
does not coincide with dk)- However, practical application of the generic formula ( |l9| ) is rather limited, because it 
does not address the full time evolution of Uk{t), and cannot strictly determine a period of modulation n/fk- 

To get an idea of how the occupation number of created particles rifc(t) evolves with time, again let us look at Fig. |l| 
and further at Fig. ^ for different values of the parameters k and q. For small and moderate q (but not too small 
k) the occupation number exhibits high frequency (period < ^) oscillations which are modulated by a long period 
behavior. For large q number of fermions jumps in a step-like manner at instances when effective mass of the fermi 
field crosses zero, superposed by very high frequency oscillations around almost constant values, as depicted in Fig. [?]. 
These jumps are modulated with a frequency ffc/2: steps in the first half of the cycle up are accumulated until nk(t) 
reaches its maximum Fk, and then steps down to zero in the second half of the cycle. 

However, this picture of high frequency features superimposed over long-period modulation is not universal. For 
small k, or for one of the most interesting cases of q 3> 1 and moderate k, the occupation number of fermions jumps 
between zero and one within time interval much shorter than period of background oscillations, as depicted by the 
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dotted curve in Fig. [1[ There are interesting situations when fermions are created in an instant (single kick) process, 
where formula (|19| ) is not applicable. Therefore, for different ranges of parameters we will shall develop different 
approaches. 



C. Semi- Analytic Theory for Averaged Occupation Number 

Numerical curves for the mode functions suggests splitting of the time evolution into higher frequency features with 
the time-scale comparable or less than the period of background oscillations < T, and low-frequency modulations 



with the period ir/vk greater than T. In this case we can utilize generic result (19). However, it is convenient to 
use not the values rik(N s T), but rather the smoothed occupation number fik(t) which is rik(t) averaged over high 

frequency oscillations, ftkir) = y J^ T+T ^ drnkir). Then we can write the smoothed occupation number of fermions 
in a factorized form 

Ufc(r) = F k sin 2 v k T , (20) 

where we introduce an envelope function Fk- The average occupation number of fermions evolves periodically with 
time. The spectrum of can be characterized by the envelope function Fk and the period of modulation 2- which 
depends also on the parameter q. 

Now we will utilize the result (Jl9|) . We found that the envelope function can be extrapolated by the factor in the 
front of sin 2 N 8 dk in ( p~S| ) and given by the expression 

F -'^fw t ( ImX ^f ■ < 21 » 

Next is to determine the frequency Vk of the fik modulations. The value dk defined after (|l^) cannot be the right 
answer, because it would incorrectly predict that the peaks of Fk are filled up first, while actually they are filled last. 
We tried the combination ///,. = tt/2 — dk, because it corresponds to correct order of saturation of Fk , and it occurs to 

work well. Therefore, the modulation frequency Vk is given by the relation cos^ fe T = -ReX^ \t). Thus, to find F k 

and Vk, one need only calculate the complex value xj^\T) after a single background oscillation, instead of performing 
a full numerical integration of Eq. jl^). 

We calculated X ( k 1] {T) numerically for ^m 2 <j> 2 background model and constructed the envelope function Fk plotted 
in Fig. ^ (similar graph was plotted in @ for ±A0 4 model). In Fig. § we show, using (^l]), how the fermionic resonance 
bands are filled after 10 background oscillations. 
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FIG. 5. The envelope function Fk for the amplitude of 
occupation number oscillations in m 2 ^ 2 -theory. Values 
for q = 10 -4 , 10 -2 and 1 are shown, from narrowest to 
broadest, respectively. 




FIG. 6. The log of the modulation period in 

m^fj) 2 -theory, where q — 1 corresponds to the middle 
curve, q = 10 -2 corresponds to the dotted curve, and 
q — 10 -4 corresponds to the remaining. 



The function i>k gives us the time scale for fermion excitation. In Fig. ^ we plot the period of modulation 2- as 
function of k. This function is peaked where Fk is peaked, i.e. the peaks of the resonance curve are the last to fill. 
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D. Method of Successive Scatterings for Fermions 



It turns out that for large values of the parameter q we can significantly advanced in calculations of nk(t) beyond 
the results of Section ( III C ). One can expect that q may be much greater that one. Indeed, q = —S 2 ■ In the context 
of chaotic inflation with the potential V(4>) — ^m 2 <f) 2 it follows from the theory of cosmological perturbations that 

,2 

S — 10 12 . On the other hand one can admit that Yukawa coupling can be h S> 10~ 6 , which provides g> 1, 

We will generalize for fermions the method of parabolic scatterings introduced for the bosonic resonance in refer- 
ence 0. The method is based on the observation that, for g> 1, change of the particle number occurs only during a 
short time interval r* near the zeros of the effective mass of the particles. This occurs because equation (^5|) for the 
number of particles in terms of the mode functions is an adiabatic invariant of the mode equation. 

For large q, the adiabaticity condition fl' k < Q 2 is violated only near the times t» when the effective mass vanishes. 
This leads to a step-like evolution in the number of fermions. This is illustrated in figure (|7|) which shows the evolution 
of nk(r) and M c g(r) = hcfjlr) + 



\ 


1 

u 

1 \ 


A 

/ \ 

! \ 

\ i ' 




i 

A i 


i 

A 

I/ \ 

J — i 




I 

I 

i - 
| 






\ i 

V 


\ , 

y 




/ 







1 2 3 4 5 

r 

2,T 

FIG. 7. nu and M c s(t) = h(f>(r) + for q = 100 and k = Jq « 31.6. 



Near the times r*, we may approximate the effective mass in equation ( |l7| ) by (m +y/qf) = (m^ + h(f>{r)) 
— 0(t*)'(t — r*) + 0((r — r*) 2 ). We may also write this as 



(m + v / 5/) « Vqf( T *)( T - T *) = Hi- ™ )( r - T *) 



(22) 



where the sign on the left hand side depends on whether /(r) is increasing (+) or decreasing (— ) at r*. Thus, in the 
neighborhood of t», the mode equation becomes 



k 2 + (q- m ) 2 (r- t*) 2 - isgn(/')y 9- m 



A fc = 



(23) 



which is a Schroedinger equation for scattering off a negative parabolic potential centered at r* . 

The method of parabolic scattering uses the exact solution of equation ( p3f ) to provide a connection formula between 
the adiabatic approximations of the full mode equation on either side of r* . Suppose the j-th zero of M c g occurs at 
time Tj. For times between Tj_i < r < Tj, the general solution of the mode equation ( |l7| ) takes the adiabatic (or 
semi-classical) form 



X{ (r) = a{N+ e- 1 So dT " fc (r) + /?£iV_ e +l £ dT < T > 



(24) 



where the coefficients cr[ and j3 k are constant for Tj—\ < t < Tj. After the scattering at Tj, Xk{t), within the interval 
Tj < t < Tj + i, again has the adiabatic form 



XI +1 (t) = al +1 N+ e- l Io drn ^ + ${+ l N^ e +l So dr ° fe(r) 



(25) 



with new coefficients a k and (3 k that are constant for Tj < t < Tj+x. At r = 0, our vacuum positive frequency 
condition requires ai = 1 and 01 = 0. Particle creation occurs when, after scattering at the times Tj, the initial positive 
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frequency wave acquires a negative frequency part. The number density of produced particles with momentum k is 
n k ~ \ flk\ 2 f° r times Tj < t < Tj+i. Furthermore, normalization requires \a k \ 2 + \f3 J k \ 2 — 1 for all j. 

The important observation is that the outgoing amplitudes (a^ +1 , can be expressed through the incoming 

amplitudes (a k , (3 k ) by means of the reflection R k and transmission D k amplitudes for scattering at tj\ 




Here 9 k = J drO(r) is the phase accumulated by the moment Tj. 
o 

If we let x — (q— m ) 1 / 4 r, equation (p3|) may be written 



(26) 



<i2 ' X " + [kl- l (-iy+x 2 ]X k = Q 1 (27) 



dx 2 

where the parameter 

2 2 

Kl ^¥Fm = ~r ==2 - (28) 

\l q-m 

A general analytic solution of equation (|27]) is a linear combination of the parabolic cylinder functions 

The reflection R k and transmission Dk amplitudes for scattering on the parabolic po- 
tential can be found from the asymptotic forms of these analytic solutions: 

Rk = ~ . = ,D k = -=== , (29) 

VI - e* A k VI - e-^l 

where the angle tfk is cpk = argT ( 1+ ^ k J + 4f U- + ln ^r^. Note that the angle if depends on the momentum k. 



Substituting (29) into (p6[), we can determine the change in induced in the a k and j3 k coefficients by a single parabolic 



scattering in terms of the parameters of the parabolic potential and the phase 9 k only. Specifically, we find 

/ Vl - e -* A le lVk -(-l)Je"t A fc+ 2 < 



j'+i 




(30) 



It is now a simple matter to find the change in particle number after one scattering. From the normalization of a k 
and Pk we have the relation 

(J -l) (^) = l«*l a -|Al a = l-2 ? 4 +1 , (31) 
which, when applied to equation ( |30| ) gives the desired result 

< +1 = e -^ + (l-2e-^,„ ; 



-2(-l) 3 e"-5 A *yi-e- 7rA 'Y / 4(l-^)sin^ ot , (32) 

where the phase 9\ ot = 29 k — ip k + arg/3^ — axgocL. Let us discuss this formula, which is the main result of our paper. 
If fermions are light, m^, — 0, their occupations number is changing with time only when background field crosses 
zeros. Without expansion of the universe, the phases 9 accumulated between successive zeros are equal. In this case 
one can try to proceed to find the solution of the matrix equation (p6j), as it was done for bosons ||. However, for 
bosons that was needed to find the stability /instability bands, which is not so interesting for fermions. In the case 
of massive fermions time intervals between successive zeros of M e ff are not equal, and problem of finding matrix 
solution of eq. ( p6| ) became even more complicated. However, in the most interesting case of expanding universe the 
phases between successive zeros of M e ff became random (see [|| for details). Then we can just put 9 to t as a random 
phase and use formula (p2|) as it is. 
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IV. PARAMETRIC EXCITATION OF FERMIONS WITH EXPANSION OF THE UNIVERSE 



To address the problem of fermionic preheating after m 2 2 chaotic inflation, we must now deal with the full 
mode equation (0), which no longer has a periodic time dependence in the complex frequency. Nevertheless, it is 
still convenient to work with the form (|l^) in the time variable r = m^t where the parameters q and k 2 are now 
understood to be time dependent: 



-2 + M cff - * 



(oMeff)' 



A(a) 



X k = 



(33) 



Check where ' stands for the derivative in respect with r, K is a comoving momentum scaled in units of m^. There 
is often used another from of the fermionic mode equation written in terms of conformal time r\ = J dt/a and mode 
function Y k = a~ 3 / 2 Xk 



dfa + [k 2 + M 2 S - id v {a M eff )] Y k = 



(34) 



This form of equation is useful, for instance, in conformal theory V(4>) = A0 4 when the problem can be reduced 
to the problem in Minkowski space-time ||. In this case background field is oscillating periodically in respect with 
time n. In the case of quadratic inflaton potential background field is oscillating with the constant period in terms 
of physical time t (or r). We therefore will use mode equation (|33|). The parameter q is now understood to be 

time dependent. Specifically, we have q{r) = 



! {r) 



scaled physical momentum p = Here, $(t) is the time 

dependent amplitude of inflaton oscillations. As is well known (c.f. Q), oscillations of <f> hi this model quickly 



approach the asymptotic solution 4>{t) 
start of inflaton oscillations, 0o ~ M 
dominated universe: a(r) = (1 + r) 2 / 3 



$(t) cos(t), $(t) 



00 



Mr, 



V37(i+t; 



where time is measured from the 



o ~ -p=, and the scale factor averaged over several oscillations behaves as in a matter 

V 07T 



A. Stochastic Parametric Excitations of Fermions 

Let us first consider the case of light fermions with a large initial resonance parameter: h<po 3> 3> ni^. Fig. (|^) 
shows a numerical solution of equation ( |i"7| ) in the absence of expansion with a resonance paramter q = 10 6 . As 
expected from for such a large q parameter, we see a step-like change in occupation number with periodic modulation. 




"A- 




FIG. 8. The occupation number n k of fermions in 
rn^ff) 2 -theory without expansion. Here q — 10 6 and the 
particular mode is k 2 ~ ^Jq. 



FIG. 9. The occupation number rik of fermions in 
rn^fj) 2 -theory with expansion taken into account. Here 
the initial resonance parameter is qo — 10 6 and the mode 



IS Kq 



In Fig. (||), we show a numerical solution to the equation ( |33] ) in the presence of expansion with the initial parameters 
qo = 10 6 . Here, the comoving number density of particles is plotted. We see that, a s for the non-expanding case, 
evolution between the zeros of the effective mass (oe equally inflaton field), M Q g w \/q{ t ) cos(t), is adiabatic. We 
can thus apply the results of Section III,D in particular the analytic formula (|32|), mutatis mutandis. The Pauli 
principle is obviously still obeyed and the typical step in particle number is still suppressed by a factor e _7rAt: 
which is now time dependent. The most important qualitative change is that the accumulated phase #tot is now 
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uncorrelated between successive zeros of the effective mass (infiaton field) . This occurs because the effective frequency 
f2ft = y/p 2 + q 2 (r) cos 2 r is no longer periodic and the accumulated phase 9 = f? + Q/, s» 2/l ^ r - > + 0(k 2 ) changes 

substantially in magnitude within one infiaton oscillation, 59^ — -£0z , after the iV s -th oscillation || . The result is that 
the sin(^tot) term in ( |32| ) becomes a random variable. As is readily apparent in Fig. (^), this destroys the periodic 
modulation of nk and the parametric excitation of fermions becomes stochastic, as anticipated in || . 

Once the periodic modulation of is destroyed, the construction of Section III,C is no longer valid: the occupation 
number cannot be characterized by an amplitude and period. In fact, the spectrum of created particles is even simpler. 
Stochastic excitation allows a given comoving mode k to obtain any amplitude in the range < n& < 1 if there are 
a sufficient number of parabolic scatterings, N s . This gives us the picture of stochastically filling a (Fermi) sphere in 
the momentum space. Numerical calculations confirm this picture, see Fig. |l^. Let us find its radius n s . A comoving 
mode will be excited if Afc(r) < 1/7T. Since Afe(r) = — -, , we have 



qo 



Therefore for the light fermions, comoving radius of excited modes is increasing with time as k s 
of the sphere is scaled as k s ~ q\j A a}^ 4 . 



(35) 



i 1 / 4 . The radius 



"A- 
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FIG. 10. The comoving occupation number of fermions in ra^4> 2 -inflation after 50 infiaton oscillations for initial resonance 



parameter q a = 10 3 . Expansion destroys the details of the resonance band and leads to a fermi-sphere of width q 1 ^ 4 which 
like a 1 / 4 while q(r) > 1. 



grows 



Due to the expansion, the amplitude $ is decreasing and the parameter q(r) drops. At some moment the Fermi 
sphere will quit expansing. Once q(r) is order of unity, the excitation is no longer strong, and redshifts of fermion 
modes will be fast enough to prevent from parametric excitation. Eventually fermions will be produced in the 
perturbative regime. If fermionic mass is non-zero, the perturbative decay goes unless 2m^ > m^. 



B. Analytic Results for Production of Supermassive Fermions 

One of the most significant differences between bosonic and fermionic parametric amplification is the possibility to 
create superheavy fermions from much liter inflatons, if inflatons are oscillating coherently. As we seen, for large q, 
massive fermions are created not when infiaton field itself crosses zero, but when the combination M e // — he/) + m^, 
crosses zero. This is because at those instances even very heavy fermions are effectively massless, because their bare 
mass is "compensated" by large value of h<f>, if <j> can have a large amplitude. This is exactly the case when (f> is the 
infiaton field. For instance, in the chaotic inflationary scenario amplitude of its oscillations immediately after inflation 
can be as large as O.lMp. 

In this case we again have a situation when mode function of heavy fermions has WKB form between zeros of M e f f , 
and can be described with parabolic scattering around this instance. Therfore our general formula ( |3^ ) works for 
superheavy fermions as well. The only modification is that the phase between successive scatterings will be defined 
by the integral over time intervals between them. In formula (t 



parameter Afc(r) will be 
A,(ri = p _ (36) 



ii 



In case of massive fermions, the criteria for excitation instead of ( J35| ) will be 

,•2 



m 



< ■ (37) 



~2 

At the start of background oscillations for large qo a term m is negligible. Thus, as in the last section, we expect 
the Fermi sphere to fill with the radius k s = ^/qoa 1 ^ 4 . However, as q drops, equation (]37]) solved for k 2 s reaches a 

/ \ 1/3 

maximum, K m . and then decreases. This maximum occurs when the scale factor is a m = I ) . This gives 

V Am J 

, x2/3 

n m — — ( q % T ) • Notice that k m scales as % . This is still not the final width of the band. This is because there 

n \4Vm/ 

are a number of background oscillations which occurs between the moment a rn and the moment when denominator in 
( |36| ) approaches zero. We put a number of oscillation before a m and between a m and the moment of zero denominator 
as approximately equal and equal to N m . In this time after a m even the modes suppressed exponentially as e _Afc 
can achieve significant occupation numbers in their random stochastic walk. We found the amplification factor in 
the from of the exponent is V 'N m . In contrary to the ligh fermions case, there is no perturbative end of the process. 

Excitation of superheavy fermions is abruptly terminated when %— m . 

We found the spectrum of created superheavy particles after the process is terminated is given by formula 

1 ( ^-l^ra) 2 \ /oca 

n k = - exp -2 - (38) 



2 "V < 

where 27 2 = In ( — ) . This formula is valid for N m greater than a few. In Fig. O we plot numerically calculated 

V 47r 2 m ) _ . — 

final spectrum of the supermassive fermions vs. an analytic formula B8, which are in a good agreement. In the case of 

— . . 1/3 

superheavy particles in an expanding universe the maximum radius of the k-sphere is scaled with qo as jk m ~ q In go, 
which is in close agreement with ||. However, one shall also check that at the moment when excitation of superheavy 
fermions terminates, their back reaction h(ipip) to the dynamics of infiaton oscillation is still negligible. 




10 20 30 40 50 60 70 



FIG. 11. Final spectrum of heavy fermions with 
masses = 50m^, when qo = 10 6 . Also shown is the 
analytic estimate (pa). 

In Fig. |ll] we plot spectrum of superheavy fermions 
with analytic formula (^). 




10 20 30 40 50 60 

K 

FIG. 12. The comoving occupation number i%k of 
fermions after M e ff crosses one zero for q = 10 6 and 
m— 50. Also shown is the analytic formula (35) which 
fits perfectly. 

t the process is terminated. Dashed curves are obtained 



C. Fermionic Production from Single Kick 

As we seen formula ( |32| ) can be extended to the case of expanding universe. Consider the very first term in formula 
(|32|). As It corresponds to the creation of fermions after effective mass M e ff(t) crosses zero for the first time. In this 
case we found 



12 



(39) 



Occupation number of fermions generated from the single kick is plotted in Fig. [l2|. Two curves corresponding to the 
formula (39) and to the numerical solutions are shown and practically indistinguishable. 



V. DISCUSSION 



In this paper we developed a non-perturbative theory for fermion production by an oscillating inflaton field. As 
we have seen, the production of fermions can be characterized as parametric excitation. Even in the simple model 
of a Yukawa coupling between fermions and a background scalar field in an expanding universe oscillating around 
the minimum of its quadratic potential, the theory of fermionic parametric excitation is rich and leads to important 
results. Fermions are created very quickly, within about ten(s) oscillations, in out-of-equilibrium states. For large 
values of the resonance parameter q = -^p, occupation number of light fermions (m^, < <C hfa) is changing 
in a step-like manner at instances tj when the inflaton amplitude (f>(tj) — passes through zero, j = 1,2,3, .. We 
have developed the method of parabolic scatterings for fermions, based closely on a similar approach for the bosonic 
resonance O. It is possible to derive a unified recursive formula, which relates the occupation number of fermions or 
bosons n J k + at the moment tj+i to the earlier value nl : 

ni +l =e-^ + (l±2e-^l lll , 



2(-l)'e-i A 2 ^l±e-^ni(l±n{)smei ot , (40) 

For bosons one shall use an upper sign and neglect the (— l)- 7 terms while for fermions one shall use the lower 
sign. For light fermions and bosons = ■ For large q the angle 9f ot can be treated as a random phase. As 

a result, formula ( p0[ ) predicts the stochastic character of parametric excitation of both bosons and fermions. In the 
fermionic case it leads to the conclusion that in the momentum space fermions chaotically fill up a broad sphere of 
the radius ~ q x ^ra. This formula also clearly shows the features of: spontaneous emission: for both bosons and 
fermions, the first oscillation leads to the spectrum nk = e _7rAfc ; stimulated emission for bosons with rik 3> 1, we have 
( n l +1 — n i) x n k'i P auu Blocking for fermions, if nl = 1, the next value will always be = 1 — e _7rAfc even in the 
stochastic case. This prevents the occupation number of fermions from exceeding 1. 

Formula (BO) can be extended to the case of massive bosons and fermions. However, here important differences 

— k 2 +m 2 

emerge. For bosons we will have A& = -7=^ , where n%b is the x-boson mass. It leads to the conclusion that 
the creation of superheavy, mj > m^, bosons is exponentially suppressed. However, for fermions we have A& = 

k 2 

- . Therefore even superheavy fermions with the mass as large as h<po can be created in abundance from 



the coherent inflaton oscillations . This occurs because the effective mass of fermions is given by the algebraic 
combination h<f>(t) + m.0 and the creation of fermions occurs when the effective mass goes to zero. The method of 
parabolic scatterings, applied for instances when this happens, leads to both formulae (^) with the corresponding 
A fc . 

There are situations where single instance (single kick) of particle creation may lead to interesting effects. An 
example is the scenario of instant preheating, which is especially important for inflationary models without minima 
of the inflaton potential fl5| . Another example is the interaction of the inflaton with superheavy fermions during 
the inflationary stage when the combination h<f>(t) + can go through zero only once. The generic formula ( |40| ) 
embraces the case where bosons and fermions are created by a single kick. We shall put j = 0, n° k = 0, and then 
n\ — e~ 7rA '' . However, instead of the parameter q we shall use another combination of coupling constants and <po and 
77i0. In this case the effect is defined by the velocity (j)* at the moment £*, which is different for bosons and fermions. 
For bosons single instance creation of particles gives the spectrum P,|l5|] 



rik = e (41) 
where t* corresponds to <j>(t*) = 0. For fermions single instance creation gives [JL4J 

n k = e h ** (42) 

where 4>(t*) + m^/h = 0. 
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We believe that the theory of fermionic preheating will be important ingredient of the realistic scenarios of the 
reheating of the universe after inflation. 
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